dat <- readstata13::read.dta13("applications/dataverse_files2/BhavnaniLeeAJPSRep1.dta")

coefFn <- function(data, w, form, coef){
  data[, 'w'] <- w
  coef(lm(form, data = data, w = w))[coef]
  
}

set.seed(02138)
form <- formula(stdlnhh100day ~ ALLaa)
coef <- 'ALLaa'
test1 <- algorithmUDP(data = dat, statistic = coefFn, B = 10, n = nrow(dat)/150, P = 150,
                     lambda = 0.35, lambda_var = 0.005, delta = 1/1300, epsilon = 1, 
                     epsilon_alpha = 1, parallelize = T, censoring_cutoff = 0.9,
                     bias_cutoff = 0.1, disjoint = T, form = form, coef = coef)


private_mod1 <- lm(form, data = dat)
private_se <- summary(private_mod1)$coef[2, 2]
private_point <- summary(private_mod1)$coef[2, 1]

plot_aa <- ggplot() + 
  geom_point(aes(y = private_point, x = 'Original'), col = 'navy') +
  geom_errorbar(aes(ymin = private_point - qnorm(0.95)*private_se, 
                    ymax = private_point + qnorm(0.95)*private_se, x = "Original"), col = 'navy', width = 0) + 
  geom_point(aes(y = test1$theta_tilde, x = 'Privatized'), col = 'dodgerblue') +
  geom_errorbar(aes(ymin = test1$theta_tilde - qnorm(0.95)*sqrt(test1$var_est), 
                    ymax = test1$theta_tilde + qnorm(0.95)*sqrt(test1$var_est), 
                    x = "Privatized"), col = 'dodgerblue', width = 0) +
  theme_bw() + 
  geom_point(aes(y = test1$theta_hat, x = 'Privatized\n(no correction)'), col = 'orange3', shape = 18, size = 4) +
 # ylim(c(0, 0.1)) + 
  scale_x_discrete(limits = c('Original', 'Privatized', 'Privatized\n(no correction)')) + 
  theme(panel.grid.major = element_blank(),  legend.position = 'none')  + 
  labs( y = 'Effect of affirmative action', x = '') + ylim(c(-0.05, 0.3))


ggsave(plot_aa, file = 'plots/fig5b.pdf', width = 3, height = 3)


dat <- readstata13::read.dta13("applications/dataverse_files2/BhavnaniLeeAJPSRep1.dta")

coefFn <- function(data, w, form, coef){
  data[, 'w'] <- w
  coef(lm(form, data = data, w = w))[coef]
  
}

set.seed(02138)
form <- formula(stdlnhh100day ~ ALLaa + as.factor(state))
coef <- 'ALLaa'
test1 <- algorithmUDP(data = dat, statistic = coefFn, B = 10, n = nrow(dat)/150, P = 150,
                      lambda = 0.5, lambda_var = 0.005, delta = 1/1300, epsilon = 1, 
                      epsilon_alpha = 1, parallelize = T, censoring_cutoff = 0.9,
                      bias_cutoff = 0.1, disjoint = T, form = form, coef = coef)


private_mod1 <- lm(form, data = dat)
private_se <- summary(private_mod1)$coef[2, 2]
private_point <- summary(private_mod1)$coef[2, 1]

plot_aa <- ggplot() + 
  geom_point(aes(y = private_point, x = 'Original'), col = 'navy') +
  geom_errorbar(aes(ymin = private_point - qnorm(0.95)*private_se, 
                    ymax = private_point + qnorm(0.95)*private_se, x = "Original"), col = 'navy', width = 0) + 
  geom_point(aes(y = test1$theta_tilde, x = 'Privatized'), col = 'dodgerblue') +
  geom_errorbar(aes(ymin = test1$theta_tilde - qnorm(0.95)*sqrt(test1$var_est), 
                    ymax = test1$theta_tilde + qnorm(0.95)*sqrt(test1$var_est), 
                    x = "Privatized"), col = 'dodgerblue', width = 0) +
  theme_bw() + 
  geom_point(aes(y = test1$theta_hat, x = 'Privatized\n(no correction)'), col = 'orange3', shape = 18, size = 4) +
  # ylim(c(0, 0.1)) + 
  scale_x_discrete(limits = c('Original', 'Privatized', 'Privatized\n(no correction)')) + 
  theme(panel.grid.major = element_blank(),  legend.position = 'none')  + 
  labs( y = 'Effect of affirmative action', x = '') #+ ylim(c(-0.05, 0.3))


ggsave(plot_aa, file = 'plots/figA4.pdf', width = 3, height = 3)